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Abstract 

We present a theoretical study of the phase diagram of a frustrated 
Ising model with nearest-neighbor ferromagnetic interactions and long-range 
(Coulombic) antiferromagnetic interactions. For nonzero frustration, long- 
range ferromagnetic order is forbidden, and the ground-state of the system 
consists of phases characterized by periodically modulated structures. At 
finite temperatures, the phase diagram is calculated within the mean-field 
approximation. Below the transition line that separates the disordered and 
the ordered phases, the frustration-temperature phase diagram displays an 
infinite number of "flowers" , each flower being made by an infinite number of 
modulated phases generated by structure combination branching processes. 
The specificities introduced by the long-range nature of the frustrating inter- 
action and the limitation of the mean-field approach are finally discussed. 

I. INTRODUCTION 

There are many physical examples in which a short-ranged tendency to order is opposed 
a long-range frustrating interactions: In diblock copolymers formed by two mutually in- 



compatible polymer chains attached to each other, the repulsive short-range forces between 
the two types of components tend to induce phase separation of the melt, but total segre- 
gation is forbidden by the covalent bonds that link the subchains together. A microphase 
separation transition occurs instead at low enough temperature, and the system then forms 
phases with a periodically modulation of structures rich in one component or the other, 
such as lamellar, hexagonal, or cubic phases Jl[ . In a similar way, self-assembly in water-oil- 
surfactant mixtures results from the competition between the short-range tendency of water 
and oil to phase separate and the stoichiometric constraints generated by the presence of 
surfactant molecules, constraints that act as the electroneutrality condition in a system of 
charged particles 0-ff]. The same kind of physics also arises in quite different fields. For 
instance, stripe formation in doped antiferromagnets like cuprates has been ascribed to a 
frustrated electronic phase separation, by which a strong local tendency of the holes to 
phase separate into a hole-rich "metallic" phase and a hole-poor antiferromagnetic phase is 
prohibited by the long-range Coulombic repulsion between the holes |5|,f|]. A last example is 
provided by the structural or topological frustration in glass-forming liquids: the dramatic 
slowing down of the relaxation that leads to the glass formation has been interpreted as 
resulting from the presence of frustration-limited domains whose formation comes from the 
inability of the locally preferred arrangement of the molecules in the liquid to tile space peri- 
odically J7|; topological frustration may also lead to low-temperature defect-ordered phases 
such as the Frank-Kasper || phases in bimetallic systems. 

A coarse-grained description of the above mentioned situations^ involves lattice or con- 
tinuum models with competing short-range and Coulombic interactions. The purpose of the 
present work is to study the phase diagram of such a model, namely the Coulomb frustrated 
Ising ferromagnet in which Ising spins placed on a 3-dimensional cubic lattice interact via 



Additional examples include cross-linked polymer mixture, interpenetrating networks, ultra-thin 
films 
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both nearest-neighbor ferromagnetic couplings and long-range Coulomb-like antiferromag- 
netic terms, the model is introduced in more detail in section [H] and its ground states as 
a function of the frustration parameter, i.e. of the ratio of the antiferromagnetic coupling 
strength over the ferromagnetic one, is studied in section |ITI[ In section [IV], we investigate 
the finite-temperature phase diagram in the mean-field approximation. Finally, the effect of 
the long-range nature of the frustrating forces, when comparing the phase behavior of the 
prototypical model with competing, but short-ranged interactions, the axial next-nearest 
neighbor Ising (ANNNI) model P-pl, as well as the limitations of the mean- field approach 



are discussed in section [V|. 

II. THE COULOMB FRUSTRATED ISING FERROMAGNET 

The model is described by the Hamiltonian 

H = -J £ SiSj + ® y £v(T ij )S i S j , (1) 

<ij> tyj 

where, J, Q > are the ferromagnetic and antiferromagnetic coupling strengths, Si = ±1 are 

Ising spin variables placed on the sites of a three-dimensional cubic lattice, < ij > denotes a 

sum restricted to nearest neighbors, r^- is the vector joining sites i and j, and v(r) represents 

a Coulomb-like interaction term with u(r)^-owhen |r| — ► oo.(In all the paper, the lattice 

spacing is taken as the unit length). 

In addition to considering the true Coulombic term, v(r) = A, we have also studied, for 

mathematical convenience in the analytical treatment, an expression of v(r) in terms of the 

lattice Green's function that satisfies the Poisson equation on the three-dimensional cubic 

lattice. In the latter case, v(r) is then simply given, up to a multiplying factor of 4tt, as the 

inverse Fourier transform of the inverse lattice Laplacian, 

, . An ^ exp(— ikr) 

V \ 2E Q=w (l-cos(A; Q ))' U 

where N is the number of lattice sites and the sum over k = (k x ,k y ,k z ) is restricted to 

the first Brillouin zone. For large |r| the lattice Green's function behaves as l/(47r|r|) so 



that the expression in Eq. (0) has the proper asymptotic behavior. In practice, even at the 
next-neighbor distance, the difference between the true Coulombic form and Eq. (Q) is very 
small fT2"fl . One has however to be careful about two points. The first one is that the v(r) 
defined in Eq. (^|) has a nonzero, finite value at r = 0, v(0) = 0.25273100986 . . . , a value 
that must of course be excluded when considering the Hamiltonian in Eq. For instance, 
Eq. ([!]) when expressed in Fourier space, can be written as 

H = ^Y,V(V\S(V\ 2 (3) 
z k 



where 

v(k) = -2 cos (^) + 4kQ 



r\- i.e.!/. i J 



v(0) 



(4) 



2J2a=x,y,z(,^ COs(k a )) 

and S(k) is the lattice Fourier transform of the Ising spin variable <S$. 

To assess the quantitative difference between the Coulombic form and that involving the 
lattice Green's function, we have calculated the energy of the system in several periodic con- 
figurations of the spins. For periodic lamellar patterns of large width m, the difference is neg- 
ligible, but it increases when m decreases. For a Neel antiferromagnetic state, the Coulombic 
energy is related to the Madelung constant and is equal to —0.873782 . . . , whereas the lattice 
Green's function expression gives —1.064356992. The second point worth mentioning is that 
the v(r) expressed in terms of the lattice Green's function has the same discrete symmetry 
as the nearest-neighbor ferromagnetic interaction whereas the true Coulombic form has con- 
tinuous rotational symmetry instead. The consequences are negligible for the ground-state 
and mean-field analyses, but may become important for the finite-temperature behavior of 



the system in other approximations |T5|. 



2 The correction term involving v(0) in Eq. (Q) has been omitted in previous papers I pjUTjj ], but 
it leads only to very small corrections when the parameter Q/J is small. 
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III. PHASE DIAGRAM AT ZERO TEMPERATURE 



In the absence of the antiferromagnetic interaction (Q = 0) the model reduces to the 
standard Ising ferromagnet, and the ground state of the system is obtained when all spins 
are aligned in a ferromagnetic state. Oppositely, when the ferromagnetic interaction is set 
to zero (J = 0), the model is equivalent to a Coulomb lattice gas and the ground-state 
of the system is a Neel antiferromagnetic state. When Q ^ 0, the Coulombic interaction 
prevents the existence of a ferromagnetic phase, and in the thermodynamic limit, the total 
magnetization (charge) is constrained to be zero. Instead, phases with modulated order, 
i.e., with periodic patterns of "up" and "down" spins subject to the constraint of zero 
magnetization, are formed. We have studied these phases both analytically with the long- 
range interactions modeled by the lattice Green's function and numerically with the true 
Coulombic form. 

For small values of the frustration parameter Q/J, the ground state consists of lamellar 
phases in which parallel planes of ferromagnetically aligned spins form a periodic structure 
along the orthogonal direction. The system in such state can mapped onto a finite one- 
dimensional system of length 2m where m denotes the width of the lamellae fl6| , |l7l . 



The short-range ferromagnetic contribution to the energy per spin of a lamellar phase 
can be readily calculated, and one gets 

E SR = -j(3--). (5) 

The Coulombic energy due to the long-range competing forces can be calculated in recip- 
rocal space by using the expression in terms of the inverse lattice Laplacian (see section [TT]). 
For a lamellar phase of period 2m the wave-vectors to be considered have only one nonzero 
component that takes the values k = ir(2n + l)/m with < n < m — 1. Correspondingly, 
the lattice Fourier transform of the Ising spin variable, S(k) = -7= exp(ikr-j), where 

N is the total number of lattice site, has its modulus given by 
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Using then the identities 

m—l j 

2 sin(7r(2n + l/(2m)) 2 = ™ ' (7) 
^ 1 m 4 2m 2 

^ sin((7r(2n + l/(2m)) 4 ~ 3 + 3 ' () 

one finds that the sum rule Ylk \S(k)\ 2 = N is properly satisfied and that the Coulombic 
energy per site, E c , is equal to 

"m—l 



E C = Q 
= Q 



. , 7T 

Ln ^o 2sin(vr(2n + l)/(2m)) 4 " 2 ™ ( } 
Tim 2 ( 1 



+ 7T - - 2u(0) 



(9) 
(10) 
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In appendix 0, we show that the same expression for the Coulombic energy is obtained 
when performing the calculation in real space with an effective one-dimensional potential, 
this latter includes a convergence factor that helps handling conditionally convergent sums 
appearing as a result of the long-range nature of the forces and that is taken to zero at the 
end of the calculation. It is worth mentioning that the above calculation in reciprocal space 
implicitly assumes that the contribution from the k = term is zero; physically, this means 
that the periodic system is embedded by a medium with infinite dielectric constant. (For 
a more detailed discussion, see Ref ||18|| .) Since for simulations on ionic systems, a similar 



choice is generally adopted fig] , both numerical and analytical calculations for our model 
have been performed using such metallic boundary conditions. 

Figure. [I] shows the total energy per spin, E/J, which is the sum of the two contributions 
in Eq. (]5|) and (|10"D for the inverse lattice Laplacian expression, as a function of Q/J for 
small values of frustration parameter Q/J. (The plot for the true Coulombic potential is 
very similar.) The slope of each straight line corresponds to the Coulombic energy and 
the intercept with the y-axis corresponds to the short-range ferromagnetic energy. For a 
given frustration parameter Q/J, the ground state is given by the straight-line that has the 
smallest energy. When Q/J^-0, the ground state consists of lamellar phases whose period 
becomes larger and larger and whose range of stability decreases. The values of {Q/J) 



corresponding to transitions between two successive ground-state structures are obtained 
by solving the equation E(m) = E(m + 1). Therefore, when Q/J goes to zero, m is 
asymptotically given by 

m ~ (Q/J)- 1/3 , (11) 

a behavior that is analogous to that observed for lamellar phases of diblock copolymer 
systems at low temperature in the so-called strong segregation limit [[0|. Note that for a 
large width m, the difference between Eq. ( |i~0l) and the value obtained fro the true Coulombic 
term is only weakly varying with m and can be taken as a constant, so that Eq. ([0]) provides 
also a very accurate estimate both for the sequence of lamellar phases and for the transition 
values of (Q/J) for the true Coulombic potential. 

When Q/J increases, the period of the lamellar phases decreases until one reaches m = 1 
for Q/J = 0.637. For a region of the frustration parameter between 0.637 and 1.800, this 
lamellar phase is then the most stable phase. In a narrow interval between Q/J = 1.800 
and Q/J = 2.122 there is a cascade of phases (1 x m 2 x oo) with m 2 decreasing until m 2 = 2 
as the frustration parameter increases (see Fig. ^|a). For 2.122 < Q/J < 3.820, the stable 
phase is (1 x 2 x oo). Note that tubular phases of type (m% x m 2 x oo) with both mi and 
m 2 > 1 are never stable. For 3.820 < Q/J < 6.237, the ground state is a 1 x 1 x oo tube. 
In another interval 6.237 < Q/J < 6.611, the system looses translational invariance in the 
third direction, and one observes a cascade of "parallepipedic" phases (1 x 1 x m 3 ), with m 3 
decreasing as the frustration increases until m 3 = 2. Between 6.611 and 9.549, the stable 
phase is (1 x 1 x 2). For Q/ J > 9.549, the ground state is the standard Neel state (lxlxl) 
(see Fig. ||a). It is noticeable that periodic structures of the type (mi x m 2 x m 3 ) with 
mi,m 2 ,m 3 > 1 and at least of the m Q 's finite always have higher energies than those of 
the sequences (mi x oo x oo), (oo x m 2 x oo), (oo x oo x m 3 ), whatever the value of the 
frustration parameter (see Appendix [BJ). 

The same analysis can be repeated with v(r) given by the true Coulomb interaction. The 
exact same sequence of ground states as before is obtained, but the values of the frustration 
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parameter at the transition points are somewhat shifted: For instance, the lamellar phase 
m — 1 is the most stable phase when 0.627 < Q/J < 5.21, the cascade of "tubular" phases 
(1 x m2 x oo) occurs around Q/J = 5.22, and (1 x 2 x oo) is stable for 5.23 < Q/J < 6.17 
(see Fig. Qb). The cascade of "orthorhombic" phases (1 x 1 x ms), with 1713 > 1, appears 
around Q/J= 14.63 (See Fig. |3|b). The standard Neel state is stable for Q/J > 15.33. It is 
worth mentioning that the counterparts of the lamellar, tubular, and orthorhombic phases in 
diblock copolymer systems [[HI (systems that are described at a coarse-grained level by the 
scalar field theory associated with the Hamiltonian in Eq. (|])) are the lamellar, columnar, 
and cubic phases, respectively. 

In the region of stability of the lamellar phases, we have also investigated if more complex 
structures involving several one- dimensional modulations of ferromagnetically ordered layers 
could be present and if multiphase transition points at which more than simply two phases 
coexist could occur. For both questions, the answer is negative. As an example, we give 
in Appendix the analytical expressions for the energy of the mixed lamellar phases that 
are formed by mixing lamellae of width m = 1 and lamellae of width m = 2 with some 
given periodic modulation. The simplest of such phases, that following the notation used 
by Fisher and Selke in their study of the axial next-nearest-neighbor Ising (ANNNI) model 
|j~0|l we denote < 1™2 P >, consists of a periodic repetition of a fundamental pattern formed 



by a succession of n lamellae of width 1 followed by p lamellae of width 2. (When n — 0, 
one recovers the simple lamellar phase of width 2, denoted < 2 >, and when p = 0, one 
recovers the simple lamellar phase of width 1, denoted < 1 >). It is easy to check that 
such mixed lamellar phases are never stable at zero temperature. In particular, they have 
a (strictly) higher energy than the pure lamellar phases at the zero-temperature transition 
point between the < 1 > and the < 2 > phases, at Q/J = 2/n, so that this latter is a simple 
two-phase coexistence point. This conclusion remains unchanged when one considers even 
more complex mixed lamellar phases, such as < \ n ^2 Pl l n2 2 P2 . . . \ ns 2 Ps >, with n a ,p a integers 
for a = 1, 2, ... s, whose fundamental period is formed by n\ 1— layer lamellae followed by 
Pi 2— layer lamellae, then by ri2 1— layer lamellae, and by P2 2— layer lamellae, etc . . . , 

8 



two successive lamellae being composed of spins of opposite signs. We have also studied 
one-dimensional quasi- crystalline arrangements of lamellae: by using binary substitution 
rules, |2(J we have built quasi-periodic structures by iteration, but we have always found 
(numerically) that their energy is higher than that of the pure lamellar phases. Note finally 
that from the calculation of the energy of the mixed lamellar phases at zero temperature, 
one can also study the change of energy induced by adding defects in the pure lamellar 



phase < 2 >. As shown in Appendix |C 2j , adding defects always cost energy, the dominant 
effect being the increase in the short-range energy. The same analysis, leading to similar 
conclusion, can be repeated for the whole range < Q/ J < 0.637 over which lamellar phases 
are favored. 

IV. MEAN-FIELD THEORY 

To describe phases with a spatial modulation of the magnetization, we consider a local 
mean-field approximation. If m; =< Si > denotes the local magnetization at site i and rh(k) 
its lattice Fourier transform, the mean-field free energy F m j is given by: 

= "2^5/ (k)l ^ (k)|2 " ^EM2cosh(^)), (12) 

where (3 = l/k B T and the effective field on site i, Hi, is equal to —JJ^jjLi Vij m ji or i n Fourier 
transformed space, 

H{k) = -JV{k)m(k) (13) 

Minimization of the free energy with respect to the local magnetizations leads to the 
self-consistent set of equations: 

rrii = tanh(PHi), (14) 

for each lattice site. One must then solve the coupled equations, Eqs. ( |T3"D and ( |T4D simulta- 
neously, and subsequently insert the solution in the expression of the free energy, Eq. (P"2"|). 
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One finally searches for the configuration of the m^'s that leads to the deepest minimum of 
the free energy for a given temperature and a given value of the frustration parameter. 

To simplify the notation in the rest of the paper, units of temperature, energy, etc . . . , 
will be chosen such that ks — J — 1- 



A. Order-disorder transition line 

Close to the transition between the disordered and the ordered phases, the magnetizations 
rrii are small and Eq. (|T4]) can be linearized, 

rrn ss (15) 

For a given value of the frustration parameter, the critical temperature T C (Q) is then given 
by 

T C (Q) = -min k y(k) = -V C (Q), (16) 

where the minimum of V^(k), V C (Q), is attained for a set of nonzero wave- vectors {k c (Q)} 
that characterize the ordering at T C (Q). For the inverse lattice Laplacian expression of the 
long-range frustrating interaction, the k c (<5)'s vary continuously with Q as follows 

k c = (± arccos(l - 0rQ), 0, 0), for < Q < 4/tt, (17) 

k c = (tt, ±arccos(3- >JnQ),0), for 4/vr < Q < 16/tt, (18) 

k c = (it, it, ± arccos(5 — \JirQ)), for 1Q/tt < Q < 36/ir, (19) 

k c = (tt, vr, tt), for 36/vr < Q. (20) 

One should of course add all vectors obtained by permuting the x, y, z coordinates in 
Eqs. (|iTD-(PU|). The above ordering wave- vectors correspond, respectively, to lamellar 
(Eq. QT7|)), tubular (Eq. QT^), orthorhombic (Eq. (0)), and cubic or Neel (Eq. phases. 
The corresponding critical temperature T C (Q) is then given by (recall that both T and Q 
are expressed in units of J) 
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T C (Q) = 6 - A^txQ + 4nQv(0) for < Q < 36/tt (21) 

T C (Q) = -6 + 4ttQ (v(0) - —) for 36/tt < Q (22) 

The mean-field approximation gives a line of second-order phase transition from the disor- 
dered to the modulated phases, with T C (Q) first decreasing with Q, reaching a minimum 
for Qmin = l/(47rt>(0) 2 ) = 1.245871 at T c (Q min ) = 6 — 1/^(0) and then increasing again for 
finally reaching a regime of linear increase with Q when Q > 3Q/ir. It is worth noting than 
the term AttQv(0) is important for large frustration: indeed, since v(0) > 1/6, it allows to 
obtain a positive critical temperature for all Q's. (In general, v(0) should be larger than 
the inverse of the number of nearest neighbors on the lattice, e.g.. 6 on a simple cubic 
lattice). For vanishing small frustrations, the critical temperature goes continuously to T c °, 
the critical temperature of the pure Ising ferromagnet. 



B. Structure combination branching processes 

At zero temperature, we showed that the system exists in pure modulated phases, whose 
modulation is commensurate with the underlying lattice. At the transition between the 
modulated and disordered phases, we have just seen that the ordering wave-vector varies 
continuously with the frustration parameter, hence indicating that a succession of incom- 
mensurate modulated phases is observed along the critical line. As the temperature increases 
from T = to T C (Q), one expects a cascade of ordered phases with commensurate spatial 
modulations of increasing complexity until a point is reached at which incommensurate 
phases appears; this is what is observed for instance in the much studied axial next-nearest 



neighbor Ising (ANNNI) model [11], in which the cascade of phases is produced by "structure 



combination branching processes" [[TTJ . To illustrate how such branching processes proceed, 
we consider first the low-temperature region of the phase diagram in which the simplest mod- 
ulated phases, the < 1 > and < 2 > lamellar phases, are stable (see Fig. |). Eqs. fll~3| ) and 
([TJ]) must now be solved beyond the linear approximation. The < 1 > phase corresponds to 
an alternate configuration of layers of up and down spins, with rrii = mi(-l) 1 , which implies 
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that the wave-vector characterizing the modulation has only one nonzero component equal 
to k = ii. Using Eq. fll^D and Eq. (|b§), one gets 

mi 



mi = tanh 
which has a nonzero solution when 



T 



[2 - ttQ + 4nQv(0)) 



(23) 



T < 2-7tQ + 47tQw(0). 



(24) 



From Eq. (|P2|), the corresponding free energies of the two phases is obtained as 



<i> 
N 



1 - ^j- + 2nQv(0) ) m\ - 7 'In 



2 cosh ( ^(2 - ttQ + AnQv{0)) 



(25) 



The < 2 > phase corresponds to an alternate sequence of pairs of ferromagnetically 
ordered layers and is characterized by a wave- vector whose only nonzero component is k = |. 
The corresponding order parameter m 2 satisfies the following equation: 



V2 



m 2 = tanh 



T 



(2 - ttQ + 2ttQv{0)) 



(26) 



which has a nonzero solution when 



T < 2-nQ + 2nQv(0). 



(27) 



The associated free energy is 
F<2> ttQ 



N 



1 



ttQv(0) ) mi -Tin 



2 cosh (^^(2 -nQ + 27iQv(0)) 



(28) 



The line of first-order transition at which the two phases coexist is defined by F < i > = F <2> . 
In the T — Q phase diagram, this line is almost vertical with Q ~ 2/V. As we have already 
stressed, no mixed phases coexist with the < 1 > and < 2 > phases at T = and Q = 2/n. 
However, the mixed < 12 > phase may become more stable at a nonzero temperature. This 
phase is characterized by the modulation m^ = ^j^-m^ cos(27ri/3 + tt/2), with the order 
parameter m\ 2 determined through the self-consistent equation 



rai2 = tanh 



^(3-^ + 4^(0) 



(29) 



12 



which has a nonzero solution for 



T < 3 



— Q + 4nQv(0). 



(30) 



The corresponding free energy is given by 



F < 



<12> 




4ttQ 



tiQv(0) m\ 



- In [2 cosh (^(3 - yQ + 4ttQt;(0)) 



AT 



3 



12 



(31) 



Comparing now the free energy of the three phases, < 1 >,< 2 >, < 12 >, one can 
show that the < 12 > phase become more stable than the other two in a wedge above 
a branching point at TJ, = 1.03 and Q = 0.63, at which the three phases coexist. This 
represents the first step of a structure combination branching process by which two adjacent 
phases, here < 1 > and < 2 >, get separated above a given branching point at finite 
temperature by a phase corresponding to the simplest combination structure, here the < 
12 > phase. A careful examination of the thermodynamic quantities shows that the entropy 
of the < 12 > phase increases more rapidly with temperature than that of the < 1 > 
and < 2 > phases; although it has an unfavorable energy contribution, the < 12 > phase 
becomes thermodynamically stable when the temperature becomes high enough (T& = 1.03). 
Moreover, for all temperatures between T = and TJ,, both the entropy and the energy of 
the < 1 > and < 2 > phases are identical along the coexistence line between the two phases. 

To study in more detail the branching processes, one can no longer rely on explicit an- 
alytical calculations as done above, because they become rapidly intractable. To obtain 
the phase diagram in a more systematic way , the mean-field equations, Eqs. flT2])- (|HD, 
can be solved iteratively for finite lattices with periodic boundary conditions. The ther- 
modynamically stable solution corresponds to that with the smallest free energy. For each 
temperature it is assumed that the spin structure repeats itself after A" layers (only com- 
mensurate lamellar structures are considered). The iterative procedure converges whenever 
the initial configuration is not too far from the equilibrium one. The iterative sequence 
is as follows. The effective fields Hi are calculated from the set of initial magnetizations 
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{m?} via Eq. (|I"4"D. One computes the Fourier transform of the fields ifj, and using Eq. ([13]) 
yields then a new set of magnetizations {m\} which is used as the input for the next it- 
eration. The calculation should be performed for various values of N for examining many 
different commensurate structures. We have carried out the calculation for N up to 16. The 
phases with large-N periodicities, like in the ANNNI phase diagram ||11|| , are only stable in 



a small neighborhood of the critical line. Therefore, the most of the phase diagram, except 
in the region close to the critical line, can be drawn by considering simple commensurate 
wave- vectors. 

The results are shown in Figs. [| and || Figure [| illustrates the structure combination 
process in the region between the < 1 > and < 2 > phases. The < 12 > and < 1 > 
phases are separated by the < 1 2 2 > phase which is itself separated at higher temperature 
from the < 1 > phase by the < 1 3 2 > phase, and so on, until one presumably reaches an 
accumulation point of the branching process along the transition line above which the < 1 > 
phase melts. Beyond this accumulation point corresponding to a sequence of < l n 2 > phases 
when n —>■ oo, devil's staircases and incommensurate phases are expected [|TTJ. Figure |5| 
provides a broader picture of the phase diagram in the region where the modulated phases 
are lamellar: the phase diagram appears as a succession of flowers of complex modulated 
phases, separated by regions in which the pure lamellar phases < 1 >, < 2 >, etc . . . , 
are stable; the flowers get closer and closer when the frustration decreases. Note that the 
first branching point at which the simplest mixed phase appears is always at a non-zero 
temperature. 

C. Soliton approach 

At high enough temperature, near the critical line, the modulated phases are incommen- 
surate since we saw in section [IV A| that the ordering wave- vector varies continuously with 



frustration. More insight in the phase behavior can be then provided by employing the soli- 
ton approach developed by Bak and coworkers |2"1~|,|2"2"|] . More precisely, this approach allows 



14 



one to study analytically the melting of a commensurate phase to incommensurate phases by 
focusing on the behavior of the domain walls that separate commensurate regions, domain 
walls that can be considered as "solitons" . In the following, we use the soliton method to 
investigate the stability of commensurate phases at high temperature. 

In the vicinity of the (upper) melting line of the < 2 > phase, one can expand the mean- 
field equations in the appropriate order variables, which leads to the following free-energy 
functional 

F 1 - 1 T 

= -7,Y,i V ( k ) + dl m k| 2 + <Hki + k 2 + k 3 + k 4 - T)m kl m k2 m k3 m k4 

T 

+ ™ 5(ki + k 2 + k 3 + k 4 + k 5 + k 6 - r)m kl m k2 m k3 m k4 m k5 ?7i k6 + . . . (32) 

dU T kl ...k 6 

where a constant term has been discarded and r is a reciprocal-lattice vector. The above 
expression, Eq. (]32|) , contains both regular and "umklapp" terms; these latter, represented 
by the second and third contributions in the right-hand side of Eq. (|32]), correspond to 
terms in which the sum of the wave- vectors is equal to a reciprocal-lattice vector, i.e., they 
keep track of the underlying lattice structure and are responsible for the stability of the 
commensurate phases. 

For studying the stability of the < 2 > phase near the critical line, we consider wave 
vectors that are close to the ordering wave vector = (tt/2, 0,0) with small fluctuations 
q x in the direction of the modulation, here along the x— axis, and q_i_ in the perpendicular 
layer. For the present case, it is sufficient to truncate Eq. (j32j) after the fourth order ||21|| , 



and after expanding the interaction term to second order in the fluctuations the free-energy 

functional can be rewritten as 

F / \ 

— = ( r + ak x + cql + c'qlj |^k, /2+q | 2 + T <5(ki + k 2 + k 3 + k 4 ) 

H=q x ,q±, qiq2Q3q4 

ri i ( 

2 mk 7 r/2+qi mk 7 r/2+q2 m -k 7r/2 +q 3 m -k^ /2 +q 4 + Y^l m k 7r/2 +qi"ik 7r/2 +q2 m k^ /2 +q3 m k 7r/2 +q4 

+ m -k V2 +qi^-k, /2 +q 2 m-k, /2 +q 3 ^-k, /2 +q 4 )] , (33) 

where r = T-A+2ttQ-AttQv(0) < 0, a = 2-2ttQ, c = 2tcQ, and d = -nQ+l. Introducing 
now two continuous order parameters, 
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m i r 



m_ r 



v/2 



d 3 q 
d 3 q 



■exp(«qr)m k +q 



^ / (2^3 exp(iqr)m_ k ^ /2+q , 



(34) 



where m_(r) is the complex conjugate of m + (r), one can express the free-energy functional 
in the following Ginzburg-Landau form: 



JV i 2 



9 . a \ 
<9x 2c / 



1 .,_ l2 1 / cr\ . l2 
+ -c\V ± m + \ + - r-- |m + | 



T T 

i I 1 4 i / 4 i 4 



(35) 



where Vj_ = (0, J^, J^) and the last term is generated by the umklapp terms. Following Bak 



and Boehm [pT| , we choose the following ansatz for the order parameters, 



m± = Aexp(±icj)(x)), 



(36) 



where the amplitude A is a constant. Note that 4>(x) is constant in the commensurate < 2 > 
phase and that A can be obtained by minimizing the free energy in that phase, which gives 
A 2 = 3\r\/T. Inserting the above expression in Eq. ( |35l) leads, up to a constant, to the 
following free-energy functional per unit area perpendicular to the x— direction: 

cA 2 



dx 



(x) <p (x) (1 — cos40(x)) 

c 12c 



(37) 



where the first term is minimized for (j)(x) = Ax /(2c), which corresponds to an incommen- 
surate spatial modulation, whereas the second term is minimized for <p(x) = 0, i.e., in the 
commensurate < 2 > phase. The overall minimization of F[(f)] is attained when the phase 
function <p(x) obeys the sine-Gordon equation: 



4>"(x) +4wsin(40(a;)) = 0, 



(38) 



where v = 7M 2 /(24c). The solution consists of regions of constant phase separated by 
solitons in which increases by | over a short distance. More generally, one can lock for a 
solution over a (large) distance L that consists of n equally spaced solitons where ^ = Lcf)' . 
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The corresponding free-energy is given by p3| , |24" 



4 i —a 



,/ 16 i, . 2n i. . . 

>' + _ V 20'exp[ — — «3 . (39) 
7T ra 



c^ 2 L 

where 0' = mr/(2L). The first term, proportional to the soliton density, is the formation en- 
ergy; the second term corresponds to a weak repulsion between solitons. The commensurate 
phase is stable as long as the first term remains positive; otherwise, the < 2 > phase becomes 
unstable with the respect to soliton formation. The critical temperature corresponding to 
this melting transition is given by 

T <2> . IC = A-2nQ + AnQv(0) - ^'J^ . (40) 

and is shown dashed line in Figs. | and |5|. 

A similar analysis can be performed to study, for instance, the stability of the < 3 > 
phase near the critical line. The main changes are that the relevant ordering wave-vector 
is now (7r/3, 0,0) and that the sixth-order umklapp terms should be kept in Eq. (|3"2"D. A 
transition between the < 3 > phase and incommensurate phases is found for 



T <3 >-/c = 5 - 4nQ + 4*Qv{0) - f II " 4 ^IV 5 f+2^Q - (41) 

The result is shown as dotted curve in Fig. [5[ The deviation from the numerical solution of 
the mean-field equations that is seen when moving further away from the critical line comes 
from the truncation of the free-energy functional that is used in the soliton approach. 

V. DISCUSSION 

The fact that spin models with competing interactions can give rise to complex spatially 
modulated phases has been known for several decades. The ANNNI model is one of the 
simplest and best-studied such system, in which Ising spin variables situated on a lattice 
are coupled via nearest-neighbor ferromagnetic interactions in (D — l)-dimensional layers 
orthogonal to, say, the and via next-nearest neighbor antiferromagnetic interactions 



along the x— axis |TTJ. The major additional ingredient that is present in the Coulomb 
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frustrated Ising ferromagnet and not is the previously studied models is the long-range na- 
ture of the competing antiferromagnetic interaction. It is then worth reviewing some of 
the differences between the phase behavior of the three-dimensional Coulomb frustrated 



ferromagnet and that of the three-dimensional ANNNI and related models ||25|| . First, the 
long 1/r range of the frustrating interaction forbids ferromagnetic ordering for any nonzero 
value of the frustration parameter Q, and, as a result, there is no Lifshitz point |26|] in the 
model. Secondly, there is no highly degenerate multiphase point at zero temperature; this 
is in contrast with the ANNNI model that possesses an infinitely degenerate multiphase 
point at zero temperature, a point from which springs an infinite number of distinct com- 



mensurate modulated phases [f27]]. Thirdly, the phase diagram of the Coulomb frustrated 
Ising ferromagnet, as illustrated in Fig. [|, displays an infinite number of distinct flowers of 
complex spatially modulated phases that emerge from zero-temperature two-phase coexis- 
tence points, the extent of the flowers in T — Q phase diagram decreasing as the frustration 
parameter decreases. 

The long range of the antiferromagnetic interaction in the Coulomb frustrated ferromag- 
net also brings about additional limitations on the mean-field approach. On general grounds, 
one can expect the mean-field description to reproduce the topology of the phase diagram 
correctly, but to become increasingly inaccurate as one approaches the transition line, both 
from below and from above, because it overlooks the role of fluctuations. In the present case, 
the fluctuations have a major effect on the transition line: as argued by Brazovskii |28] 



on 



the basis of a self-consistent Hartree treatment of a field-theoretical model analogous to the 



Coulomb frustrated Ising ferromagnet, and confirmed by Monte Carlo simulations |13|, the 
fluctuations drive the order-disorder transition from second to first-order. The mean-field 
approach is thus questionable in the vicinity of the transition line (which is why it may 
not be worth pursuing the search for devil's staircases as was done for the ANNNI model 



2l|,|29[|). However, the main points reviewed above are not affected. 
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FIG. 1. Ground-state energy of the lamellar phases with m = 1,2, ..7 as a function of the 
frustration parameter Q/J for the inverse lattice Laplacian expression of v(r). When Q/J goes to 
zero, the period of the lamellar phases increases. The inset zooms in on the region of vanishing 
frustration parameters which corresponds to the left down region of the figure. 
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1 2 3 4 51234567 

Q/J Q/J 

FIG. 2. Ground-state energy of the tubular phases with m = 1, 2, ... 7 (full lines) and m = oo 
(dashed line) as a function of the frustration parameter Q/J. (a) For the inverse lattice Laplacian 
expression, the tubular phases, (1 x mi x oo) with m<i > 1„ are stable for frustration parame- 
ters 1.800 < Q < 3.820 (b) For the true Coulombic potential, the tubular phases are stable for 
5.21 < Q < 6.17. < 1 >, < 1 x 2> and < 1 x 1 > are short-hand notations for < 1 x oo x oo >, 
< 1 x 2 x oo > and < 1 x 1 x oo >, respectively. The vertical dot-dashed lines are visual guidances 
for denoting the stability region of phases. 
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FIG. 3. Energy of the orthorhombic phases with m = 1,2, ... 7 as a function of the frustra- 
tion parameter Q/J. (a) For the inverse lattice Laplacian expression, the orthorhombic phases 
(1 x m<i x oo) with ni2 > 1, are stable for frustration parameters 6.237 < Q < 9.549 (b) For the 
true Coulombic potential, the orthorhombic phases are for 14.63 < Q/J = 15.33. < 1 x 1 > is a 
short-hand notations for < 1 x 1 x oo >. The vertical dot-dashed lines are visual guidances for 
denoting the stability region of phases. 
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FIG. 4. Mean-field phase diagram: structure combination branching processes occurring at 
finite temperatures for a region of frustration parameter where the ground states consists of < 1 > 
and < 2 > phases. The dashed line corresponds to the soliton-approach prediction for the (upper) 
stability of the < 2 > phase. The units are chosen such that ks = J = 1. 



24 




0.2 0.4 0.6 0.8 1 

Q 

FIG. 5. Mean-field phase diagram: a (partial) view of the infinite sequence of "flowers" of com- 
plex modulated phases appearing at finite temperatures for the range of the frustration parameter 
where the ground states are lamellar phases. The dashed and dotted lines correspond to the soli- 
ton-approach prediction for the (upper) stability of the < 2 > and < 3 > phases, respectively. The 
units are chosen such that ks = J = 1. 
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APPENDIX A: CALCULATION OF THE COULOMBIC ENERGY OF 
LAMELLAR PHASES IN REAL SPACE 



The calculation of the ground-state energy due to the long-range Coulombic interaction 
can be performed as follows. For lamellar phases, the sum over the reciprocal vectors is per- 
formed along one direction. One then obtains the one- dimensional potential corresponding 
to the inverse lattice Laplacian [ |T2]j , 

W{i) = 2nQ eM ~ a ^-\ (Al) 
a 

where a is a convergence factor that will taken to zero at the end of the calculation. The 
introduction of a > in the calculation allows one to handle conditionally convergent sums. 
The average Coulombic energy per site of a cell 2m is given by 

i +oo 2m 

E c = 7- E Y i W{\i-j + 2mN'\)S i S j , (A2) 

4m N'=-oo i,j 

where N' is the index for labelling the right and left cells and i,j are denote sites within 
the cell. The energy per site can be divided into two parts: the first one comes from the 
interaction between the spins within the cell and the second one comes from the interaction 
between a spin and its images in the other cells; this reads 

Eq = Ei + E 2 

-i 2m i 2m 

= 7-E^(|i-« + r £ £W(|z-j + 2miV')^, (A3) 



By considering all contributions between pairs of sites within the cell, one gets for E 



1 



m— 1 2m— 1 



Ex = m(W(0) - W{m)) + J2 ( 2 ( m - n) - n)W{n) - (( 2m ~ n ) ~ n)W{n) . 

2m \ n=l n=m+l J 

(A4) 

After some calculation, and taking the limit a — > at the end, one gets 

_ / 2irm 2 7r\ . , . 

E ^ Q — + 3 ■ (A5) 
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The sum over the right and the left cells can also be performed, and E2 is then given by 

ttO „-a(2m+i-j) , -a(2m+j~i) 

E 2 = ^ty] i SiSj. (A6) 

2m a(l - e - a2m ) 3 y ' 

By using the electro-neutrality condition (J2i Si = 0), and in the limit a — > 0, one obtains 
i? 2 ) which has a finite value: 

„r\ 2m 

1 v^.-et 12 



2m 2 fr{ 
TiQm 2 



(A7) 



This gives for the Coulombic energy E c 



inn 2 7r 



E < = Q — + l ■ (A8) 



By subtracting the self-energy of the inverse lattice Laplacian potential to Eq. ( |A8|) , one 
exactly recovers Eq. fllPl). 

APPENDIX B: GROUND-STATE ENERGY OF TUBULAR AND 
ORTHORHOMBIC PHASES 

Let us calculate the energy per site for configurations whose phases are periodic with a 
orthorhombic cell (mi x m 2 x 777.3). The short-range energy per spin is obtained as 

/ 2 2 2 \ 

E SR = -J(3 (Bl) 

V mi m 2 m^J 

In the reciprocal space, the allowed wave-vectors have components ((2ni + l)iv/mx, {^ n 2 + 
l)7r/m2, (2n 3 + l)7r/m 3 ), with < n\ < mi — 1, < n% < m 2 — 1, < n 3 < m 3 — 1, so that 
the Coulombic energy per spin E c is given by 



E C = Q 



mi — 1 ni2 — 1 m:j — 1 . 

, o /(in _liW\2, o / to*, _i_iw\2 27117(0) 



-=o -=o "3=0 \(ELi sin ) ITLi sin ( 



(2n a +l)7r 
2m, 



(B2) 
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When both m 2 and go to infinity, one obtains lamellar phases of period 2m, and 
Eq. ([BSD educes to Eq. ([10]). 

When the periodic structure looses translational invariance in a second direction, one 
obtains a lattice of tubes whose perpendicular section is a rectangular cell m\ x 1712. Although 
we have not obtained a fully analytical expression for such phases, some results can be 
derived. The Coulombic energy for tubes of section mi x m2 is bounded as follows: 

E c (inf(mi,m 2 )) < E c (m 1 ,m 2 ) < E c (sup(m 1 ,m 2 )). (B3) 

If p x mi = m 2 = m where p is a positive integer, the energy E c (m) behaves as 

E c (m 2 ) ~ C(p)m 2 + 0(m), (B4) 

where the numerical coefficients C(p) are summarized in table Q. 

We have also calculated numerically, via Eqs. ( pl[) and (|B2|) the total energy of modulated 
orthorhombic phases (m x x m 2 x m 3 ) when 00 > m 1; m 2 , m 3 > 1. It always higher than that 
of the phases (1 x m 2 x 777.3) or (1 x 1 x 7713). 

APPENDIX C: MIXED LAMELLAR PHASES 

The inverse Laplacian approximation allows to calculate exactly the energy of a large 
number of periodic structures at zero temperature. As an illustration we present here the 
results for mixed lamellar phases that are potential ground-state candidates in the region 
of Q/J where the most stable among the simple lamellar phases involve lamellae of width 
m = 1 (phase < 1 >) and m = 2 (phase < 2 >). 

1. The < l n 2 p > phases. 

The < l n 2 p > phases are the simplest mixed phases that one can construct with the two 
elementary bricks formed by lamellae of width m = 1 and m = 2: they are formed by a 
periodic sequence of ferromagnetically aligned layers whose fundamental period consists of 
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n one-layer lamellae followed by p two-layer lamellae, two successive lamellae being formed 
by spins of opposite signs. Because of the electoneutrality (zero magnetization) condition 
one has to distinguish three different families: (n = 2q, p = 2r), (n — 2q + 1, p — 2r) and 
(n = 2q, p = 2r + 1), where q and r are integers. The < l 2 <i+ 1 2 2r+1 > phase is not allowed at 
zero temperature because they do not satisfy the electoneutrality condition. Because of the 
two-dimensional in-layer ferromagnetic ordering, the wave-vectors characterizing lamellar 
phases have only one nonzero component. 



a. The < l 2q 2 2r > phases. 



The size of the one-dimensional unit cell is L = 2q + 4r and the allowed values of 
the nonzero components of the wave vector are k = 2 ^ l 4r where / is an integer such that 
/ = 0, . . . , 2q + 4r — 1. Summing over all sites of the lattice, one finds 



S(k) 





sin(2rA;) 




2q + 4r 


cos(/c)cos(| 


)l 



N 



sin(g/c) 



' k\ I ' 



if k 7^ ±| and k ^ n. Otherwise, one obtains | S(±^) \— 2 ^ r v / 8r and | S(tt) 



2q+4r 



(CI) 



2q. 



Using the identities 



g+2r- 1 

E 



and 



1=1 tDb \ q +2r) t0b \2(q+2r)) 



q +2r-l S in(^-) 2 

^ sinf^P =r(r + g); 

1=1 bU1 \ q+ 2r) 



(C2) 



(C3) 



one can express the Coulombic energy per spin as 

2q 2 + 16r 2 + 16gr 



E C = Q 



7T- 



2ttv(0) 



(C4) 



(2q + 4r) 2 

The short-range contribution can be easily calculated and the total energy per spin is equal 
to 

2g 2 + 16r 2 + 16gr 



E 



<1 2q 2 2r > 



J -2 + 



g 

g + 2r 



Q 7T- 



(2? + 4r) J 



2ttv(0) . 



(C5) 
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It is now easy to show that the above energy, whatever the strictly positive values of q and r 
and whatever the value of Q/J, cannot be less than either the energy of the < 1 > phase or 
that of the < 2 > phase. Indeed, for the conditions E <1 2 q2 2r > < E < \ > and E <1 2 q2 2r > < E <2> 
to be simultaneously satisfied, one must have 

< Q/J < 4^T, (06) 

ixq n(q + r) 

which is impossible. 

b. The < l n 2 2m > and < l 2n 2 m > phases. 

To satisfy the requirement of the global electroneutrality, the unit cell of such 

phases is built as follows. For the < i™2 2m > phases, the unit cell is 

Tltl • • • T IITT • • • IITT III • • • I TTII • • • TTII an d the nonzero components of the wave- 
„ . . „ ' 

n m n m 

vector are given by k = (g^+^w) wnere I — 1 ■ ■ ■ 8m + 2n — 1. After some algebra, the Fourier 
transform of the spin variable S(k) is obtained as 



m i= ',t: 1 2 A , (ct) 





2 


cos(f ) 




2n + 8m 


cos(fc)cos(|) 



if k 7^ 7r and | S(tt) \— 2 //+ Sm 2n. ^ ne total energy per spin is then 

f 2n \ „ ( 2n 2 + QAm 2 + 32mn ,.\ . , 

The < Y 2n 2 m > phase is characterized by the se- 

quence tltl • • • TITTII • • • IITT IT • • • IT II • • • TTTII an d the nonzero components of the 

V v *• v ' S v ' S v ' 

n m n m 

wave- vector are given by k = j^^) w ith 1=1... 4m + 4n — 1. This leads to 

I S(k) I- ^ 2 1 Sin(nk) 1 (C9) 
4n + 4m|cos(A;)co S (|) | 1 J 

if k 7^ ±|, and | S(±^ ) |= A ^^ n 2\[2n. The total energy per spin is finally given by 

n + mj ' ~* y (4n + 4m) 2 

(CIO) 

As before, there is no range of the frustration parameter Q/J for which these phases become 
ground states of the system. 

30 



n \ „ ( 32n 2 + 64m 2 + 32mn n , A 
E <ln2 2 m> = J -2 + + Q\ n — 2nv{0) . 



2. Energy of defects in the < 2 > phase. 



The issue of the stability of the mixed versus simple lamellar phases can be addressed 
in a different way. Starting from the < 2 > phase, one can calculate the energy change 
brought by inserting one or several defects of type < 1 > within the periodic structure. The 
creation of one such defect results from the flip of a pair of up-down spins. This corresponds 
to the simplest excitation that one expects in the < 2 > phase: ... TTIITTII ••• =^ 
• • • TTITIT-LI • • • The resulting defective structure can also be viewed as a < l 4 2 2m > phase 
where m — > oo. The Coulombic energy of a < l 4 2 2m > structure for a cell of size 4 + 4m is 
derived from Eq. (|C4|) and reads 



e c< i4 2 2ra> = (4 + Am)E c<1 4 2 2 m> = Q 



2 + 4m 2 + 8m . 

7T 2tto(0)(4 + 4m) 

1 + m v yv ' 



(Cll) 

If one subtracts the Coulombic energy of the < 2 > phase for the same unit cell (/ = 4 + 4m) 
from the above equation, one obtains the Coulombic energy for one excitation in a < 2 > 
phase. Strikingly, this energy goes to zero when m goes to infinity, which means that the 
presence of one defect in the < 2 > phase does not change at all the Coulombic energy of 
this phase. (This is not a result of the macroscopic limit.) Conversely, this defect has a 
short-range energy cost that is easily obtained as Esr = 4. In all the region of frustration 
parameters Q/J where the < 2 > phase is more stable than any other simple lamellar phase, 
the presence of one defect is then energetically unfavorable for the system. 

We have also calculated the energy of a < 2 > phase where two defects have been intro- 
duced. We have first calculated the energy of the < 2 2m l 2p 2 2n l 2 " > and < 2 2m+2n l 2p+2q > 
phases which are identical. Setting p = q = 2 whereas n,m — > oo, one obtains the energy 
of the < 2 >phase in the presence of two defects. We have also found that the introduction 
of the two defects is energetically unfavorable for the system for any value of the frustra- 
tion. Although we have not obtained a general proof, the phase-locking into simple lamellar 
phases at zero temperature seems to be well established for the model. 
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TABLES 

TABLE I. Coefficient of the leading term of the Coulombic energy per site for periodic tubular 
structures p = 1, 2, 3, 5, oo 



P Cip) C( P )*p 2 

1 0.22208 0.22208 

2 0.089083 0.35633 

3 0.045957 0.41361 

4 0.027568 0.44109 

5 0.018304 0.4576 

oo vr/6 = 0.52359... 
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